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ABSTRACT 

In the context of count-in-cells statistics, the joint probability distribution of the 
density in two concentric spherical shells is predicted from first first principle for 
sigmas of the order of one. The agreement with simulation is found to be excellent. This 
statistics allows us to deduce the conditional one dimensional probability distribution 
function of the slope within under dense (resp. overdense) regions, or of the density 
for positive or negative slopes. The former conditional distribution is likely to be more 
robust in constraining the cosmological parameters as the underlying dynamics is less 
evolved in such regions. A fiducial dark energy experiment is implemented on such 
counts derived from ACDM simulations. 
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1 INTRODUCTION 

With the advent of large galaxy surveys (e.g. SDSS and in 
the coming years Euclid (iLaureiis et al.ll201lh . LSST), as¬ 
tronomers have ventured into the era of statistical cosmol¬ 
ogy and big data. Hence, there is a dire need for them to 
build tools that can efficiently extract as much information 
as possible from these huge data sets at high and low red- 
shift. In particular, this means being able to probe the non¬ 
linear regime of structure formation. The most commonly 
used tools to extract statistical information from the ob- 
serv ed galaxy distribut i on ar e N-point correlation functions 
(e.g IScoccimarro et aH 1 19981 ) which quantify how galaxies 
are clustered. In our initially Gaussian Universe the matter 
density field is fully described by its two-point correlation 
function. However departure from Gaussianity occurs when 
the growth of structure becomes non-linear (at later times or 
smaller scales), providing information that is not captured 
by the two-point correlation function but is recorded in part 
in the three-point correlation function. Obviously N-point 
correlation functions are increasingly difficult to measure 
when N increases. They are noisy, subject to cosmic vari¬ 
ance and highly sensitive to systematics such as the complex 
geometry of surveys. It is thus essential to find alternative 
estimators to extract information from the non-linear regime 
of structure formation in order to complement these classical 
probes. This is in particular critical if we are to understand 
the origin of dark energy, which accounts for ~ 70 % of the 
energy budget of our Universe. 

One such method to accurately probe the non-linear 


regime is to implement perturbation theory in a highly 
symmetric configuration (spherical or cylindrical symmetry) 
for which the full joint cumulant generating functions can 
be constructed. Such constructions take advantage of the 
fact that non-linear solutions to the gravitational dynamical 
equations (the so-called spherical collapse model) are known 
exactly. Corresponding observables, such as galaxy counts in 
concentric spheres or discs, then yield very accurate analyt¬ 
ical predictions in the mildly non-linear regime, well beyond 
what is usually achievable using other estimators. The corre¬ 
sponding symmetry implies that the most likely dynamical 
evolution (amongst all possible mapping between the initial 
and final density field) is that corresponding to the spherical 
collapse for which we can write an explicit linear to nonlinear 
mapping. This has been demonstrated in the l imit of zero 
variance u sing direct diagram resummations llBernardeaul 
ll99lll994lB which was later shown to correspond to a sad¬ 
dle approximation lljuszkiewicz et al.l 1 19931 : IValageaj l2002l : 
iBernardeau et al.ll20l4l . The key point on which this whole 
paper is based upon, is that the zero variance limit is shown 
to pro vide a remarkably good working model for finite vari¬ 
ances llBernardeaull 1994 : IBernardeau et al.ll20l3 ). 

This formalism also allows to weigh non-uniformly dif¬ 
ferent regions of the universe making possible to take into 
account the fact that the noise structure in surveys is not 


^ The original derivations were actually derived from the hierar- 
chical model that aimed a t describing the fully nonlinear regime, 
llBalian &: SchaeffeilllQSsI) 
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homogenous. For instance, low density regions are probed 
by fewer galaxies. Conversely, on dynamical grounds, we 
also expect the level of non-linearity in the field to be in- 
homogenous: low density regions are less non-linear. Hence 
it is of interest to build statistical estimators which probe 
the mildly non linear regime and that can be tuned to probe 
subsets of the field, offering the best compromise between 
these constraints. In the context of the cosmic density field, 
the construction of conditional distributions naturally leads 
to the elaboration of joint probability distribution functions 
(PDF hereaft er) of the density in conc entric cells. 

Following iBernardeau et all l|2014h (hereafter BPC), we 
propose in Section 2 to extend one-point statistics of den¬ 
sity prohles and to the full joint probability distribution 
function of the density in two concentric spheres of differ¬ 
ent radii. This is obtained using perturbation theory core 
results on the cumulant generating function, the double in¬ 
verse Laplace transform of which is then computed from 
brute force numerical integration. From that PDF, we will 
also present the statistics of density profiles restricted to 
underdense (resp. overdense) regions, and the statistics of 
density restricted to positive (resp. negative) slopes (Sec¬ 
tion 3). Theoretical predictions will be shown to be in very 
good agreement with simulations in the mildly non-linear 
regime. Dependence with redshift will also be discussed. Fi¬ 
nally Section 4 presents a simple hducial dark energy exper¬ 
iment, while 5 wraps up. 


2 THE 2-CELL DENSITY STATISTICS 

For the sake of clarity, let us present and briefly comment the 
formalism. We consider two spheres Si of radius Ri {i = 1,2) 
centered on a given location of space xq. Our goal is to derive 
the joint PDF of the density in iSi and 1 S 2 denoted pi and 
rescaled so that {pi) = 1. 


2.1 The cumulant generating function 


In the cases we are interested in, the joint statistical prop¬ 
erties of pi and p 2 are fully encoded in their moment gener¬ 
ating function 

°° 

Mr^rF>^i,\2) = ^ (P?P2)^-y, (1) 

= (exp(Aipi -I- A 2 P 2 )) , (2) 


that can be related to the cumulant generating 
function, (Ai, A 2 ), through At(Ai, A 2 ) = 

exp ^ 2 )], so that 

exp [<pHifl 2 ('^ii -^ 2 )] 

= j dpidp2pBifl2(Pi.P2)exp(Aipi-I-A 2 P 2 ), (3) 

where Vr^rFF^ P 2 ) is the joint PDF of having density pi in 
5i and p 2 in 52. We will now exploit a theoretical construc¬ 
tion that permits the explicit calculation of 'Pflii ?2 (pi, P 2 ). 


2.1.1 Upshot 

As we will sketch in the following, this theoretical construc¬ 
tion yields the explicit time dependence of the Legendre 


transform of PhiHs (Ai, A 2 ) in the quasi-linear regime. Such 
a Legendre transform is defined as 

^'fllH 2 (Pl,P 2 ) = Aipi -I- A 2 P 2 — ¥^flifl 2 (Al, A 2 ), (4) 


where p; are determined implicitly by the stationary condi¬ 
tions 

3 

A = ^4'HiH2(pi,P2), 1 = 1,2. (5) 

The fundamental relation is then that, in the limit of zero 
variance, this Legendre transforms taken at two different 
times, 4t(pi,p2;?7) and 4''(pi, p 2 ; p^), take the same value 

^'flifl2(Pi>P2;p) = (p'i,P2;p') , (6) 

provided that piRf = p'iR'^, and that p' and pi are linked 
together through the nonlinear dynamics of spherical col¬ 
lapse. 

Equation m, when applied to an arbitrarily early time 
Tj', yields a relation between 'L(pi,p 2 ;p) and the statistical 
properties of the initial density fluctuations. In particular, 
for Gaussian initial conditions, ^{pi, p 2 \rji) can easily be 
calculated and expressed in terms of elements of covariance 
matrices, 

«'flifl2(Pi>P2;Pi) = ^ X] 

ijsS2 


where Hij is the inverse of the matrix of covariances, Eij = 
(TiTj), between the initial density contrasts in the two con¬ 
centric spheres of radii Ri. One can then write the cumulant 
generating function at any time through the spherical col¬ 
lapse mapping between one final density at time p in a sphere 
of radius Ri and one initial contrast in a sphere centered on 
the same point and with radius R) = Rip\^^ (so as to en¬ 
compass the same total mass); it can be written formally 
as 


pi = Csc{v,n) 


1 

(1 - D+{ri)T/n)‘' ’ 


( 8 ) 


where, for the sake of simplicity, we use here a simple pre¬ 
scription, with D+{'q) the linear growth factor and n = 
21/13 to reproduce the high-z skewness. 

Recall that only (pi,P 2 ) is easily computed. 

The statistically relevant cumulant generating function, 
<Pflifl 2 ('^ii -^ 2 ), is only accessible via equation (|4]) through an 
inverse Legendre transform which brings its own complica¬ 
tions. In particular note that all values of Ai are not accessi¬ 
ble due to the fact that the pi - Ai relation cannot always be 
inverted. This is signaled by the fact that the determinant of 
the transformation vanishes, e.g. Det [9pi9pj'I'(pi, P 2 )] = 0. 
This condition is met both for finite values of pi and Ai. The 
corresponding contour lines of ip(Ai, A 2 ) was investigated in 
BPC and successfully compared to simulation. 


2.1.2 Motivation 

It is beyond the scope of this letter to re-derive equa¬ 
tions dl-® - a s omew hat deta iled presentation c an be 
found in IValagea^ (1200 2l ) and in iBernardeau et ahl (l2014h 
- but we can give a hint of where it comes from: it is always 
possible to express any ensemble average in terms of the sta¬ 
tistical properties of the initial density held so that we can 
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formally write 

exp [(/?]= Jt>tiT>T 2 V{ti, T 2 ) exp (Aipi(ri) + X 2 P 2 ir 2 )) ■ (9) 

As the present-time densities pi can arise from different ini¬ 
tial contrasts, the above-written integration is therefore a 
path integral (over all the possible paths from initial condi¬ 
tions to present-time configuration) with measure 'Dt\'Dt 2 
and known initial statistics P(ri,r 2 ). Let us assume here 
that the initial PDF is Gaussian so that, 

f>(n.n)d„d„ = '^°Pl-*<".-")l dnd„ , (10) 

Itt 

with 'F then a quadratic form. 

In the regime where the variance of the density field is 
small, equation m is dominated by the path correspond¬ 
ing to the most likely configurations. As the constraint is 
spherically symmetric, this most likely path should also re¬ 
spect spherical symmetry. It is therefore bound to obey the 
spherical collapse dynamics. Within this regime equation 
becomes 

exp [(/3]~ jdTidT 2 P(ri, r 2 ) exp (AiCsc(n) + A 2 Csc('r 2 ))(ll) 

where the most likely path, pi = Csc(’ 7 ,'ri) is the one-to- 
one spherical collapse mapping between one final density 
at time rj and one initial density contrast as already de¬ 
scribed. The integration on the r.h.s. of equation m can 
now be carried by using a steepest descent method, ap¬ 
proximating the integral as its most likely value, where 
Aipi(ri) -I- \ 2 P 2 {t 2 ) — ^{ti,T 2 ) is stationary. It eventually 
leads to the fundamental relation dSl) when its right hand 
side is computed at initial time (and the fact that ([5]) is 
valid for any times p and r]' is obtained when the same rea¬ 
soning is applied twice, for the two different times). 

The purpose of this letter is to confront numerically 
computations of the two-cell PDF derived from the expres¬ 
sion of <p(Ai,A 2 ) with measurements in numerical simula¬ 
tions. 


2.2 The 2-cell PDF using inverse Laplace 
transform 


Once the cumulant generating function is known in equa¬ 
tion ®, the 2-cell PDF, 'P{pi,p 2 ), is obtained by a 2D in¬ 
verse Laplace transform of <p(Ai, A 2 ) 



dAi 

27ri 


L 


^ exp(-^ Pi\i -(- (p(Ai, A2)), 
i^2^i 


( 12 ) 


with ip given by equations @-®. From this equation, it 
is straightforward to deduce the joint PDF, V{p, s), for the 
density, p = pi and the slope s = {p 2 — pi)Ri/AR, AR 
being R 2 — Ri, as 

^/i +‘fiiKu)), (13) 


with A = Ai -f A 2 , p = X 2 AR/Ri. Following this dehnition, 
(p(A,p) is also the Legendre transform of ^{pi,s — {p 2 — 
pi)Ri/AR). 

In order to numerically compute equation (HD , we sim¬ 
ply choose the imaginary path (Ai,A 2 ) = i(niAA,n 2 AA) 
where ni and n 2 are (positive or negative) integers and the 


2-cell PDF 

6 . 
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p = density 

Figure 1. Joint PDF of the slope (s) and the density (p) as 
given by equation inj for two concentric spheres of radii Ri = 10 
Mpc/h and R 2 = 11 Mpc/h at redshift 2 : = 0.97. Dashed contours 
corresponds to Log R = 0, —1/2, —1, • • • — 3 for the theory. The 
corresponding measurements are shown as a solid line. 



s=slope 

Figure 2. Dependence of the PDF of the slope on the number 
of points used in the numerical integration in m- The reference 
PDF is computed using 1000^ points (dark blue) and is compared 
to the result of the numerical integration when using 320^ (blue) 
and 480^ (light blue) points. 


step AA has been set to 0.15. The maximum value of Xi used 
here is 75 resulting into a discretisation of the integrand 
on 1000^ points. Fig. [T] compares the result of the numer¬ 
ical integration of equation m to simulations. The corre- 
sponding dark matter simulation (carried out with Gadget2 
(ISpringei[200il l is characterized by the following ACDM 
cosmology: Dm = 0.265, Da = 0.735, n = 0.958, Ho = 70 
km-s“^-Mpc“^ and as = 0.8, fib = 0.045 within one stan - 
dard deviation of WMAP7 results llKomatsu et al.ll201ll L 
The box size is 500 Mpc/h sampled with 1024® particles, 
the softening length 24 kpc/h. Initial conditions are gener- 
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Figure 3. Density profiles in underdense (solid light blue), over- 
dense (dashed purple) and all regions (dashed blue) for cells of 
radii Ri = 10 Mpc/h and i ?2 = 11 Mpc/h at redshift z = 0.97. 
Predictions are successfully compared to measurements in simu¬ 
lations (points with error bars). 


Figure 4. Density PDF in negative slope (solid light blue), pos¬ 
itive slope (dashed purple) and all regions (dashed blue) for cells 
of radii Ri = 10 Mpc/h and R 2 = 11 Mpc/h at redshift z = 0.97. 
Predictions are successfully compared to measurements in simu¬ 
lations (points with error bars). 


ated using mpgraf ic llPrunet et al.|[200^ . An Octree is built 
to count efficiently all particles within concentric spheres of 
radii between ii = 10 and llMpc//i. The center of these 
spheres is sampled regularly on a grid of 10Mpc//i aside, 
leading to 117649 estimates of the density per snapshot. 
Note that the cells overlap for radii larger than 10 Mpc/h. 

The convergence of our numerical scheme is investigated 
by varying the number of points. Fig. [5] shows that the nu¬ 
merical integration of the slope PDF has reached 1% preci¬ 
sion for the displayed range of slopes. Obviously, the integra¬ 
tion is very precise for low values of the slope and requires 
a largest number of points for the large-slope tails. 


3 CONDITIONAL DISTRIBUTIONS 

3.1 Slope in sub regions 

Once the full 2-cell PDF is known, it is straightforward to 
derive predictions for density profiles restricted to under¬ 
dense 


■P{s\p < 1 ) = 


Jo dp PjFs) 

/.Zo^s/o dp:P(p,s) ’ 


and overdense regions 


P(s|p > 1) 


JJ° dp r{p,s) 
fZo ds fo dp V{p, s) 


(14) 


(15) 


Fig.|3]displays these predicted density profiles in underdense 
and overdense regions compared to the measurements in our 
simulation. A very good agreement is found with some slight 
departures in the large slope tail of the distribution. As 
expected, the underdense slope PDF peaks towards posi¬ 
tive slope, while the overdense PDF peaks towards negative 
slope. The constrained negative tails are more sensitive to 
the underlying constraint, providing improved leverage for 
measuring the underlying cosmological parameters. 


3.2 Density in regions of given slope 


Conversely, one can study the statistics of the density given 
constraints on the slope. For instance, the density PDF in 
regions of negative slope reads 


P{p\s < 0) 


/o°° dpjJ^dsV{p,s) 


(16) 


Fig. |4] displays the predicted density PDF in regions of pos¬ 
itive or negative slope. As expected, the density is higher 
in regions of negative slope. An excellent agreement with 
simulations is found. 


3.3 Redshift evolution 

Fig. [ 5 ] displays the density profiles in underdense and over- 
dense regions as measured in the simulation for a range of 
redshifts. This figure shows that the high density subset for 
moderately negative slopes is particularly sensitive to red¬ 
shift evolution, which suggests that dark energy investiga¬ 
tions should focus on such range of slopes and regions. 


4 FIDUCIAL DARK ENERGY EXPERIMENT 

Let us conduct the following fiducial experiment. Con¬ 
sider a set of 10,000 concentric spheres, and measure 
for each pair the slope and the density, {pi,Si}. Recall 
that the cosmology is encoded in the parametrization of 
the spherical collapse on the one hand (v), and on the 
linear power-spectrum, Pj’^, (via the covariance matrix, 
= J Pj'"{k)W{Rik)W{Rjk)<Pk/{2-Kf with W{k) = 
3 (sin(fc)/fc — cos(fc)) /k^) on the other hand. For scale in¬ 
variant power spectra with power index n, given equa¬ 
tion (15]l , we have a three parameter (n, u, a) set of mod¬ 
els. For a parametrized PDF, Pn,v,<j{p,s) given by equa¬ 
tion m, we can compute the log likelihood of the set as 
C{n,v,a) = dog Pn,^,a{pi, Si). Fig. in] displays the corre¬ 
sponding likelihood contours at one, three and five sigmas in 
the simple case in which only one parameter {a here) varies. 
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s=slope 

Figure 5. Same as Fig. |3] for a range of redshifts as labeled. 
Only the underdense (p < 1) and the overdense (p > 1) PDFs are 
shown. 



Figure 6. Log-likelihood for a fiducial experiment involving 
10,000 concentric spheres of 10 and 11 Mpc//i measured in our 
simulation. The model here only depends on the variance cr [v 
and n are fixed). The contours at 1, 3 and 5 sigmas centered on 
the true value 0.23 are displayed with dark blue dashed lines. The 
same experiment can be carried out when the three parameters 
vary. 

This experiment mimics the precision expected from a sur¬ 
vey of useful volume of about (350h“^Mpc)® which is found 
to be at the percent level. This work improves the findings 
of BPC which relies on a low-density approximation for the 
joint PDF. 


5 CONCLUSION 

Extending the analysis of BPC, predictions for the joint 
PDF of the density within two concentric spheres was 
straightforwardly implemented for a given cosmology as a 
function of redshift in the mildly non-linear regime. The 
agreement with measurements in simulation was shown on 
Figs. [T] |3] and |4] to be very good, including in the quasi- 
linear regime where standard perturbation theory normally 
fails. A fiducial dark energy experiment was implemented on 
counts derived from ACDM simulations. 

Such statistics will prove useful in upcoming surveys as they 


allow us to probe differentially the slope of the density in re¬ 
gions of low or high density. It can serve as a statistical indi¬ 
cator to test gravity and dark energy models and/or probe 
key cosmological parameters in carefully chosen subsets of 
surveys. The theory of count-in-cells could be applied to 2D 
cosmic shear maps so as to predict the statistics of projected 
density profiles. Velocity profiles and combined probes in¬ 
volving the density and velocity fields should also be within 
reach of this formalism. 
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